Geostatistical integration of rock physics, seismic amplitudes and geological models in North-Sea turbidite systems

Turbidite systems are complex heterogeneous siliciclastic deposits calling for a new approach to geophysical and geostatistical reservoir characterization. The aim of this paper is to demonstrate the integration of seismic amplitude and rock physics data into realistic geological scenarios using a novel geostatistical approach. Avseth (2000, 2001 with others) presented an integrated rock physics/seismic inversion approach to predict seismic facies within the Glitne Field reservoir, North Sea. The outcome of this work is a seismic derived conditional probability of facies. We provide a geostatistal methodology for integrating this large-scale probability model with smaller scale geological models consisting of channel complexes typical for known turbidite systems around the world. Introduction Turbidite reservoirs currently represent major hydrocarbon targets in several areas in the world. These deep-water clastic systems are often characterized by a complex sand distribution, stretching the limits of current, conventional seismic and geostatistical modeling and analysis tools. Nevertheless, reservoirs of this type must produce at high rate in order to return the large drilling and production cost. Hence the reservoir heterogeneity must be accurately quantified and the associated uncertainty measured in order to determine the capital investment risk involved. We propose a methodology to create fine-scale reservoir models of turbidite systems, constrained by pre-stack seismic and well-log data. We then apply it to a field in the preproduction phase in the North-Sea. In the first stage, a probability model of observing sandy and shaly facies – as originally presented by Avseth et al. (2001, Geophysics) and Avseth (2000, PhD Dissertation) – is determined based on the seismic data only. The original contribution of this paper comes in the second phase, when this coarse-scale seismic derived probability map is integrated with smaller scale variations of submarine channels using a new geostatistical method. The second step is needed because 1) the deep seismic (2km) cannot detect individual smaller scale channels which may act as conduits of fluids, 2) the seismic model is not locally constrained by small scale well-log data and 3) although the seismic derived facies probablity model provides uncertainty on the absence or presence of facies, it does not provide uncertainty of future cumulative oil production. The geostatistical approach proceeds first by simulating a reservoir training image using a Boolean simulation algorithm for channels. This training image is not constrained by any reservoir specific data, it is merely conceptual. Next a pixel-based geostatistical simulation method uses this training image to constrain a set of alternative facies models (gridblocks of size 12.5x12.5x1m) to the seismic derived facies probability, the small scale well data and the geometric patterns of channels as depicted by the training image. We show that our methodology is fast, general and integrates consistently all the available data at relevant scales. Conditional distributions of facies from seismic The Glitne Field To showcase our approach, we use data from the Glitne field, a reservoir of Tertiary age, located in the South Viking graben in the North Sea. A comprehensive database is available with well logs from 7 wells, including P-wave velocity (Vp), shear wave velocity (Vs), density, gamma ray and resistivity data. All wells have Vp logs, only two wells have Vs logs. Pre-stack amplitude data from 2D lines and a 3D seismic amplitude cube is available. Facies associations in turbidite systems In deep-water siliciclastic systems one can encounter six different lithofacies at the seismic scale. Table 1 provides a definition for each facies type. However in the Glitne field only 4 facies types, namely type II, III, IV and V are recognized, representing a gradual transition from clean sandstone to pure shale. Moreover, facies II is subdivided into three classes, namely IIa, IIb and IIc, see Table 1. These facies can be linked to depositional sub-enviroments and sedimentary processes within the deep-water clastic system (Walker, 1978). An idealized sketch of a submarine fan system is depicted in Figure 1. The upper fan is characterized by channel fill turbidite conglomorate systems, debris flow or slump deposits (facies I and VI) or by starved shale units (facies V). The turbidite currents on the upper fan are usually transported through a single deep channel depositing conglomorated and thick bedded sands (facies I and II). This feeder-channel is usually confined by overbanks. In the mid-fan and lower-fan areas coarse grained material is transported radially via (lobe)channels (facies II), while more fine grained and interbedded deposits are typical in inter-channel and marginal areas of the lobe (facies III and IV). Figure 1 represents a very idealized, conceptual model and a greater complexity is expected in most systems due to the influence of different controlling factors (e.g. dominant grain-size, type of the feeder system, tectonic activity, paleobathymetri). Modelling facies probability distributions from seismic In this section, the method for modeling the conditional distribution of facies, from seismic data developed in Avseth (2000) is summarized. The seismic amplitude data A(u), is available on a 3D grid of N nodal locations u, u 2 D. The goal is to model the distribution of each of the K facies types represented by an indicator variable I(u; sk), k = 1; : : : ;K where, I(u; sk) = 1 if class sk occurs 0 else We define the following classes s1 = facies IIa s2 = facies IIb s3 = facies IIc s4 = facies III s5 = facies IV s6 = facies V The target is to model the conditional distribution of facies at each location ui PrfI(ui; sk) = 1 j fA(u);8u 2 Dgg; i = 1; : : : ; N; 8sk (1) Note that the conditioning is not only on the co-located seismic amplitude datum A(u), but on all seismic data (or at least a restricted set of amplitudes). Amplitude data is a result of reflection on contrasting geological events, hence the information of facies at ui contained in amplitude data is spread laterally and vertically (Yilmaz, 1987). Avseth (2000) propose to determine the conditional probability (1) in five steps, summarized as follows: 1. Using well-log data from a single type well, one models the conditional distribution of facies given co-located Vp and gamma-ray data. 2. The conditional probability in 1. is used to perform a classification of facies at other wells. The classification result is used to determine for each facies, the distribution of seismic properties such as impedance and Vp=Vs ratio. This provides the variability of seismic properties per facies. 3. The conditional distributions of seismic properties obtained in 2. are corrected for fluid presence with the BiotGassmann theory. 4. The distribution of seismic properties per facies, allows quantifying the joint distribution of zero-offset reflectivity Ro versus AVO gradient G for each of the facies classes. 5. The seismic amplitude data A allows calculating Ro(u) and G(u) at each location u. Hence the joint distribution calculated in 4. can be used to determine the target conditional probability of (1). We elaborate each step in more detail. 1. The Vp and information from a single type-well are used to calibrate a facies classification method. Using a detailed geological interpretation of the well-log and core data, the facies succession in the type well are determined within the framework of Table 1. This provides a training data set of facies at each depth location and the well-log Vp and at the same depth locations. A non-parametric density classification method is used for this purpose. The success rate of the classification model is higher than 80 %. Figure 2 shows crossplots of density versus gamma ray and Vp versus gamma ray from the type-well. 2. Using the classification system calibrated from the single type-well, the vertical facies succession in other wells is determined on the basis of their Vp and logs. Performing this classification at all remaining 6 wells provides us with a probability distributions of seismic properties, i.e. acoustic impedance and Vp=Vs ratio, for each of the six facies types. This distribution is denoted as PrfZj(u) z j i(u; sk) = 1g k = 1; : : : ;K j = 1; 2 (2) where Z1(u) is the acoustic impedance and Z2(u) the Vp=Vs ratio. 3. Seismic properties are also dependent on the particular fluids present in each block of the reservoir. Only two fluid types are considered: oil and brine. The fluid type is represented by an indicator variable j(u; tk) with t1 is oil and t2 for brine. Using the Biot-Gassmann theory, the conditional distributions (2) for the sandy facies (facies II and III) are corrected for oil based on the observed brine saturated values. We denote these conditional distributions as PrfZj(u) z j i(u; sk) = 1; j(u; tk0) = 1g; (3) k = 1; : : : ;K k0 = 1; 2. 4. Given the variability of seismic properties per facies, one can establish the expected variability of seismic amplitude data as follows. First, the number of facies and fluids classes are redefined: c1 represents oil-sands (oil+facies IIa,b,c, III), c2 represents brine-sands (brine+facies IIa,b,c, III), c3 represents shales (oil+brine+facies IV, V). Facies IV (category s5) is a typical cap-rock of the reservoir under study. Therefore, half-space models are created with facies IV in the upper half and each of the other facies classes c1, c2 and c3 in the lower half. These half-space models allow to analytically calculate the Amplitude-Versus-Offset response (AVO). In general, AVO seeks to extract rock parameters by analyzing seismic amplitudes variations as function of the reflection angle. For angles less than 30 degrees incidence, one finds that the reflectivity for half-space models equals R(u; ) = R(u; 0) +G(u) sin( ) with R(u; 0) termed the zero-offset reflectivity and G(u) the reflectivity gradient, both depending on the seismic properties R(u; 0) = 1 2 Vp(u) Vp(u) + (u) (u) G(u) = 1 2 Vp(u) Vp(u) 2 V 2 s (u) V 2 p (u) 2 Vs(u) Vs(u) + (u) (u) V (u) is the difference between the velocity at u and the velocity of the facies one block above u, in this case the facies IV. the difference in density. Using the distribition of seismic properties per facies, a Monte Carlo simulation allows determining the joint variation of R(u; 0) and G(u) for each of the half-space models, which allows estimating PrfI(u; ck) = 1 j R(u; 0); G(u)g k = 1; : : : ; 3 (4) 5. Finally, the pre-stack seismic amplitude information of the top layer of the reservoir (the layer under the cap-rock with facies IV) allows to determine R(u; 0) and G(u) for each grid block in that top layer. Using a least-square procedure, the parameters R and G are extracted from the amplitude data. Since one knows the probability (4) one can determine, for given R(u; 0) and G(u) calculated from the seismic amplitudes, the conditional distribution (1) for each of the three facies types. These conditional probabilities are two-dimensional. A 3D extension of this probability model has also been constructed (Mukerji et al., 2001), but not used in this study. Figure 3 shows the seismic derived probability for the oil sand facies lumped together with an the interpreted feeder channel and turbidite fan system. Next we show how fine-scale geological models can be combined with this seismic derived probability model to construct reservoir models for flow simulation and reservoir management purposes. Geostatistical integration of seismic and geological model Seismic provides a large scale model of the presence of sand facies. The outcome of the previous section is a conditional probability of facies in blocks whose lateral resolution is in the order of 100m (extend of the Fresnel zone at these depth). However the actual CDP sampling (not to be confused with resolution) results in 12:5 12:5m blocks. The vertical resolution of the 2D probability model is around 15m The resulting seismic derived conditional probability can be used to perform classification of the facies in each block individually, however the resulting such classification model has clear limitations The conditional probability (1) provides only a singlepoint measure of uncertainty, i.e. it measure the uncertainty of each grid-block one at a time, independent from other grid-blocks. It is an estimated facies model that might not reflect the finer scale geological understanding of the deposit. Although the R and G parameters are calibrated from fine scale well data, the resulting conditional probability (1) does not necessarly honor the well data at the finer scale. Well data are used as a basis for obtaining the AVO pdfs, but any fine scale well-log data is missing in the seismic derived probability. In this section, we use a novel geostatistical methodology to integrate the seismic derived conditional probability of facies with a given geological model. We concentrate on the lowerchannel/upper-fan area indicated in Figure 3. This particular part of the reservoir is not penetrated by any well. Geological models Submarine channel systems have been documented from many ancient successions world-wide, the present-day seafloor and in ongoing hydrocarbon exploration and production. Unlike fluvial channel systems, limited quantitative information is available due to limited outcrop data at the reservoir modeling scale, often strong compaction and structural deformation and the bias towards detecting only the largest channels (Clark and Pickering, 1996). The scale of seismic data is usually larger than 15 m in the vertical, hence quantitative measurements of smaller scale features (1m vertical scale) are almost non-existant in the literature. Nevertheless, the purpose of this report is to provide a novel methodology for integrating any arbitrary geological model with large scale seismic data. In fact the proposed method allows to deal with the potentially large uncertainty in the geological models The single feeder channel system in Figure 3 consists of various smaller scale architectural features. In outcrop, one often observes a stacking of channel sand bodies, demonstrating channel growth resulting from the interaction between lateral and vertical amalgamation processes (Figure 4). The latter being caused by changing basin floor topograhy, salt or muddiapir action or synsedimentary faulting. Evidently, the channel growth has a strong control on the interconnectivity and the width-to-depth ratio of reservoir sand bodies. Figure 4 shows a diagram adapted from Clark and Pickering (1996) of various stacking patterns numbered 1 to 9. In mud-rich systems the development of levees will often promote a strong degree of vertical aggradation, whereas in sand-rich systems channels can expect to have a strong degree of lateral migration due to the absense of the mud-rich levees. The Glitne field has been classified in class 6 in Figure 4, which has a good vertical connectivity, a somewhat lesser lateral connectivity. From statistical data (Clark and Pickering, 1996; Weimer et al. 1995), it can be concluded that single channels with the feeder channel complex, have a width to thickness ratio from roughly 1:10 to 1:50, and most have a low to medium sinuousity (See Clark and Pickering (1996) for a definition) of 1 (straight) to 1.5. Individual channels can be anywhere between 50m to 250m wide, while the entire feeder-channel complex can be more than 1km wide. In a first approach, we will model only two facies: channel and non-channel, where the sand facies is represented by the indicator I(u) = 1 if sand occurs in gridblock u, equals zero else. This means that the oil-sand and brine-sand are combined into a single facies category. In a second phase (not in this work), one could model in more detail the within-channel variation, by including more detailed architectural elements such as levees. Methodology First, we use an unconditional Boolean simulation method to construct a so-called analog reservoir. The analog reservoir, also termed training image, depicts the believed variability of channels within the lower-channel/upper-fan part of the deposit. The training image is not conditioned to any seismic or well data, it is a purely conceptual model. Next, we discretize this single boolean model on a pixelgrid. The reservoir analog is used to extract the channel patterns and reproduce them in a set of simulation models conditioned to the seismic derived probability (1) and well-data (if available). This approach is markedly different from applying a direct conditional Boolean simulation. A direct Boolean simulation would have difficulty handling the seismic derived probability of facies (1) which is defined on a pixel-grid. Even the more advanced Boolean methods (Viseur, 2000) can only handle vertical and areal proportion information for each facies. Interpreting the seismic derived conditional probability (1) as an areal proportion would be a too limited interpretation of the seismic data. Moreover, the Boolean approach is not general, can handle only objects and is difficult to condition to multiple well data (not applicable in this particular case), and requires considerable CPU. The proposed pixel-based approach does not suffer from either of these factors. The training image approach allows the geologist to interactively adjust the facies architecture in the training image, a concept which is appealing in the modeling of complex systems such as turbidites. Multiple training images could be used, each reflecting different scales of variability or complexity of the heterogeneity. Training image need not carry any reservoir specific information, which makes them easier to construct; they reflect a prior geological concept. In this study, we limit ourselves to a single large training image. We use the fluvsim algorithm (Deutsch and Wang, 1996) to construct unconditional Boolean simulations. The proposed approach therefore takes the best of both world: it relies on the geological realism of the Boolean methods, yet retains the flexibilty and speed of the pixel-based approaches. From training image to conditional distribution Traditional to geostatistics is to capture the geological variability in a variogram. A variogram is measure for the degree of correlation/connectivity between any two single locations in space. Since the variogram is only a two-point statistics, it cannot model curvi-linear structures such as channels. The representations of complex geological features requires the modeling of so-called multiple-point statistics. Two approaches exists to the modeling of these multiple-point statistics (Caers, 2000, Strebelle and Journel, 2000, 2001). We use the Strebelle and Journel-approach which is faster in 3D. The idea of the socalled multiple-point geostatistics is to infer measures of correlation between multiple spatial locations (multiple-point statistics), from the training image or reservoir analog. The Strebelle and Journel approach is termed snesim (single normal equation simulation) and classifies under the sequential simulations algorithms in geostatistics. The general outline of the sequential simulation method is as follows Construct a fine 3D grid of N nodes and assign well-data to closest grid nodes Define a random path Until each non-datum node ui on the random path is visited step 1. Search for closest nearby data and previously simulated nodes. step 2. Construct a conditional distribution for each facies based on the searched set of nodes. step 3. Draw from the conditional distribution a facies category and assign that simulated category to node ui. end simulation Traditional to geostatistics is to use kriging in step 2., considering each datum one at a time to construct the conditional distribution. However, kriging relies on the variogram, hence traditional geostatistical methods cannot reproduce complex geological models. In the snesim method one models directly the conditional distribution from the training image, no kriging or variogram modeling is involved. The conditioning data, consisting of well data and previously simulated nodes are pooled together into one single data event. The conditional distribution in step 2. is obtained by scanning the training image for replicates of that particular single data event. In our case, only two facies categories are present, hence only the sand-conditional probability is needed. This conditional distribution is simply the proportion of times a sand facies (or shale) occurs in the center location of the data event. This approach would require scanning a potentially large training image each time a node u is to be simulated, since the geometrical configuration changes from node to node. The snesim algorithm solves this problem by pre-scanning (before the until loop starts) for all combinations of data events (n + i 1) and storing them in a dynamically allocated search tree (see Strebelle and Journel, 2000, for details). Data integration We can extend the snesim method easily to allow for integration of the seismic derived conditional probability (1) with the geological model. First a frequentist approach is given to the probability (1) as follows. If L reservoirs, denoted as i(u);u 2 D, l = 1; : : : ; L are simulated using the snesim algorithm, we should expect that in the extension of the snesim approach the proportion of sand pL(u) at each node calculated over the L realizations approximates the conditional probability (1).